#region Using Directives

using System;
using System.IO;
using NUnit.Framework;
using dnA.Math;
using dnA.Exceptions;

#endregion Using Directives

namespace dnA.dnAtests.Math
{
	
	// suite of tests for complex double symmetric algorithm
	[TestFixture]
	public class ComplexDoubleSymmetricLevinsonTest
	{

		// unit testing - case 1
		
		ComplexDoubleMatrix T1;
		ComplexDoubleMatrix Inverse1;
		ComplexDouble Det1;
		ComplexDoubleVector X1;
		ComplexDoubleVector Y1;
		ComplexDoubleVector A1;
		
		// unit testing - case 2
		
		ComplexDoubleMatrix T2;
		ComplexDoubleMatrix Inverse2;
		ComplexDouble Det2;
		ComplexDoubleVector X2;
		ComplexDoubleVector Y2;
		ComplexDoubleVector A2;
		
		// unit testing - case 32
		
		ComplexDoubleMatrix T3;
		ComplexDoubleMatrix Inverse3;
		ComplexDouble Det3;
		ComplexDoubleVector X3;
		ComplexDoubleVector Y3;
		ComplexDoubleVector A3;
		
		#region Text Fixture Setup

		[TestFixtureSetUp]
		public void SetupTestCases()
		{
			
			//            +-                           -+
			//            |  2 + 2 I,    1,      1/2    |
			//            |                             |
			// T1 =       |     1,    2 + 2 I,    1     |
			//            |                             |
			//            |    1/2,      1,    2 + 2 I  |
			//            +-                           -+
			
			T1 = new ComplexDoubleMatrix(3);
			T1[0, 0] = T1[1, 1] = T1[2, 2] = new ComplexDouble(2.0, 2.0);
			T1[1, 2] = T1[1, 0] = T1[2, 1] = T1[0, 1] = new ComplexDouble(1.0, 0.0);
			T1[2, 0] = T1[0, 2] = new ComplexDouble(0.5, 0.0);
			
			// Det1 =     - 39/2 + 23/2 I

			Det1 = new ComplexDouble(-39.0/2.0, 23.0/2.0);
			
			//            +-                                                                 -+
			//            |  223/1025 - 289/1025 I,   1/82 + 9/82 I,   - 23/1025 + 39/1025 I  |
			//            |                                                                   |
			// Inverse1 = |      1/82 + 9/82 I,     31/164 - 49/164 I,     1/82 + 9/82 I      |
			//            |                                                                   |
			//            |  - 23/1025 + 39/1025 I,   1/82 + 9/82 I,   223/1025 - 289/1025 I  |
			//            +-                                                                 -+
			
			Inverse1 = new ComplexDoubleMatrix(3);
			Inverse1[0, 0] = Inverse1[2, 2] = new ComplexDouble(223.0/1025.0, -289.0/1025.0);
			Inverse1[1, 1] = new ComplexDouble(31.0/164.0, -49.0/164.0);
			Inverse1[1, 2] = Inverse1[1, 0] = Inverse1[2, 1] = Inverse1[0, 1] = new ComplexDouble(1.0/82.0, 9.0/82.0);
			Inverse1[2, 0] = Inverse1[0, 2] = new ComplexDouble(-23.0/1025.0, 39.0/1025.0);

			//            +-   -+
			//            |  1  |
			//            |     |
			// X1 =       |  2  |
			//            |     |
			//            |  3  |
			//            +-   -+

			X1 = new ComplexDoubleVector(3);
			X1[0] = new ComplexDouble(1.0, 0.0);
			X1[1] = new ComplexDouble(2.0, 0.0);
			X1[2] = new ComplexDouble(3.0, 0.0);
			
			//            +-            -+
			//            |  11/2 + 2 I  |
			//            |              |
			// Y1 =       |    8 + 4 I   |
			//            |              |
			//            |  17/2 + 6 I  |
			//            +-            -+
			
			Y1 = new ComplexDoubleVector(3);
			Y1[0] = new ComplexDouble(11.0/2.0, 2.0);
			Y1[1] = new ComplexDouble(8.0, 4.0);
			Y1[2] = new ComplexDouble(17.0/2.0, 6.0);
			
            //               +-                  -+
            //               | -29/130 + 14/65 I  |
            // A1 =          |                    |
            //               |   -8/65 + 1/65 I   |
            //               +-                  -+

			A1 = new ComplexDoubleVector(2);
			A1[0] = new ComplexDouble(-29.0/130.0, 14.0/65.0);
			A1[1] = new ComplexDouble(-8.0/65.0, 1.0/65.0);

			//            +-                                                            -+
			//            |        4,         12 - 4/3 I,  80/3 + 64/3 I,  48 - 16/3 I   |
			//            |                                                              |
			// T2 =       |    12 - 4/3 I,        4,         12 - 4/3 I,  80/3 + 64/3 I  |
			//            |                                                              |
			//            |  80/3 + 64/3 I,   12 - 4/3 I,        4,         12 - 4/3 I   |
			//            |                                                              |
			//            |   48 - 16/3 I,  80/3 + 64/3 I,   12 - 4/3 I,        4        |
			//            +-                                                            -+
			
			T2 = new ComplexDoubleMatrix(4);
			T2[0, 0] = T2[1, 1] = T2[2, 2] = T2[3, 3] = new ComplexDouble(4.0, 0.0);
			T2[1, 0] = T2[2, 1] = T2[3, 2] = T2[0, 1] = T2[1, 2] = T2[2, 3] = new ComplexDouble(12.0, -4.0/3.0);
			T2[2, 0] = T2[3, 1] = T2[0, 2] = T2[1, 3] = new ComplexDouble(80.0/3.0, 64.0/3.0);
			T2[3, 0] = T2[0, 3] = new ComplexDouble(48.0, -16.0/3.0);
			
			// Det2 =     - 119677696/81 - 71916032/81 I
			Det2 = new ComplexDouble(-119677696.0/81.0, -71916032.0/81.0);
			
			// Inverse2 = array(1..4, 1..4,
			//            (1, 1) = - 4078146891/1189860020660 + 1198583631/594930010330 I,
			//            (1, 2) = 841024329/237972004132 + 647975751/237972004132 I,
			//            (1, 3) = 681375750/59493001033 - 2744980005/118986002066 I,
			//            (1, 4) = 1627068543/594930010330 + 4330930839/594930010330 I,
			//            (2, 1) = 841024329/237972004132 + 647975751/237972004132 I,
			//            (2, 2) = - 4050233403/237972004132 + 1620441/59493001033 I,
			//            (2, 3) = 3651059103/237972004132 + 7203556509/237972004132 I,
			//            (2, 4) = 681375750/59493001033 - 2744980005/118986002066 I,
			//            (3, 1) = 681375750/59493001033 - 2744980005/118986002066 I,
			//            (3, 2) = 3651059103/237972004132 + 7203556509/237972004132 I,
			//            (3, 3) = - 4050233403/237972004132 + 1620441/59493001033 I,
			//            (3, 4) = 841024329/237972004132 + 647975751/237972004132 I,
			//            (4, 1) = 1627068543/594930010330 + 4330930839/594930010330 I,
			//            (4, 2) = 681375750/59493001033 - 2744980005/118986002066 I,
			//            (4, 3) = 841024329/237972004132 + 647975751/237972004132 I,
			//            (4, 4) = - 4078146891/1189860020660 + 1198583631/594930010330 I
			//            )
			
			Inverse2 = new ComplexDoubleMatrix(4);
			Inverse2[0, 0] = new ComplexDouble(-4078146891.0/1189860020660.0, 1198583631.0/594930010330.0);
			Inverse2[0, 1] = new ComplexDouble(841024329.0/237972004132.0, 647975751.0/237972004132.0);
			Inverse2[0, 2] = new ComplexDouble(681375750.0/59493001033.0, -2744980005.0/118986002066.0);
			Inverse2[0, 3] = new ComplexDouble(1627068543.0/594930010330.0, 4330930839.0/594930010330.0);
			Inverse2[1, 0] = Inverse2[0, 1];
			Inverse2[1, 1] = new ComplexDouble(-4050233403.0/237972004132.0, 1620441.0/59493001033.0);
			Inverse2[1, 2] = new ComplexDouble(3651059103.0/237972004132.0, 7203556509.0/237972004132.0);
			Inverse2[1, 3] = Inverse2[0, 2];
			Inverse2[2, 0] = Inverse2[0, 2];
			Inverse2[2, 1] = Inverse2[1, 2];
			Inverse2[2, 2] = Inverse2[1, 1];
			Inverse2[2, 3] = Inverse2[0, 1];
			Inverse2[3, 0] = Inverse2[0, 3];
			Inverse2[3, 1] = Inverse2[0, 2];
			Inverse2[3, 2] = Inverse2[0, 1];
			Inverse2[3, 3] = Inverse2[0, 0];

			//            +-   -+
			//            |  4  |
			//            |     |
			// X2 =       |  3  |
			//            |     |
			//            |  2  |
			//            |     |
			//            |  1  |
			//            +-   -+
			
			X2 = new ComplexDoubleVector(4);
			X2[0] = new ComplexDouble(4.0, 0.0);
			X2[1] = new ComplexDouble(3.0, 0.0);
			X2[2] = new ComplexDouble(2.0, 0.0);
			X2[3] = new ComplexDouble(1.0, 0.0);

			//            +-                 -+
			//            |  460/3 + 100/3 I  |
			//            |                   |
			// Y2 =       |   332/3 + 40/3 I  |
			//            |                   |
			//            |    488/3 + 80 I   |
			//            |                   |
			//            |     300 + 40 I    |
			//            +-                 -+
			
			Y2 = new ComplexDoubleVector(4);
			Y2[0] = new ComplexDouble(460.0/3.0, 100.0/3.0);
			Y2[1] = new ComplexDouble(332.0/3.0, 40.0/3.0);
			Y2[2] = new ComplexDouble(488.0/3.0, 80.0);
			Y2[3] = new ComplexDouble(300.0, 40.0);
			
			//            +-                                    -+
			//            |   -500663/1194095 - 177563/170585 I  |
			//            |                                      |
			// A2 =       |       -11884/2191 + 7762/2191 I      |
			//            |                                      |
			//            |  + 399838/1194095 - 328742/170585 I  |
			//            +-                                    -+

			A2 = new ComplexDoubleVector(3);
			A2[0] = new ComplexDouble(-500663.0/1194095.0, -177563.0/170585.0);
			A2[1] = new ComplexDouble(-11884.0/2191.0, 7762.0/2191.0);
			A2[2] = new ComplexDouble(399838.0/1194095.0, -328742.0/170585.0);

			//            +-                         -+
			//            |   5,  4 I,  3,  2 I,  1   |
			//            |                           |
			//            |  4 I,  5,  4 I,  3,  2 I  |
			//            |                           |
			// T3 =       |   3,  4 I,  5,  4 I,  3   |
			//            |                           |
			//            |  2 I,  3,  4 I,  5,  4 I  |
			//            |                           |
			//            |   1,  2 I,  3,  4 I,  5   |
			//            +-                         -+

			T3 = new ComplexDoubleMatrix(5);
			T3[0, 0] = T3[1, 1] = T3[2, 2] = T3[3, 3] = T3[4, 4] = new ComplexDouble(5.0, 0.0);
			T3[1, 0] = T3[2, 1] = T3[3, 2] = T3[4, 3] = T3[0, 1] = T3[1, 2] = T3[2, 3] = T3[3, 4] = new ComplexDouble(0.0, 4.0);
			T3[2, 0] = T3[3, 1] = T3[4, 2] = T3[0, 2] = T3[1, 3] = T3[2, 4] = new ComplexDouble(3.0, 0.0);
			T3[3, 0] = T3[4, 1] = T3[0, 3] = T3[1, 4] = new ComplexDouble(0.0, 2.0);
			T3[4, 0] = T3[0, 4] = new ComplexDouble(1.0, 0.0);
			
			// Det3 =     2160
			
			Det3 = new ComplexDouble(2160.0, 0.0);

			//            +-                                             -+
			//            |   17/60,  -1/10 I,  -4/15,   1/15 I,   7/60   |
			//            |                                               |
			//            |  -1/10 I,   1/5,   -1/30 I,  -2/15,   1/15 I  |
			//            |                                               |
			// Inverse3 = |   -4/15,  -1/30 I,   7/15,  -1/30 I,  -4/15   |
			//            |                                               |
			//            |   1/15 I,  -2/15,  -1/30 I,   1/5,   -1/10 I  |
			//            |                                               |
			//            |    7/60,   1/15 I,  -4/15,  -1/10 I,  17/60   |
			//            +-                                             -+

			Inverse3 = new ComplexDoubleMatrix(5);
			Inverse3[4, 4] = Inverse3[0, 0] = new ComplexDouble(17.0/60.0, 0.0);
			Inverse3[4, 3] = Inverse3[3, 4] = Inverse3[1, 0] = Inverse3[0, 1] = new ComplexDouble(0.0, -1.0/10.0);
			Inverse3[4, 2] = Inverse3[2, 4] = Inverse3[2, 0] = Inverse3[0, 2] = new ComplexDouble(-4.0/15.0, 0.0);
			Inverse3[4, 1] = Inverse3[1, 4] = Inverse3[3, 0] = Inverse3[0, 3] = new ComplexDouble(0.0, 1.0/15.0);
			Inverse3[4, 0] = Inverse3[0, 4] = new ComplexDouble(7.0/60.0, 0.0);
			Inverse3[3, 3] = Inverse3[1, 1] = new ComplexDouble(1.0/5.0, 0.0);
			Inverse3[3, 2] = Inverse3[2, 3] = Inverse3[2, 1] = Inverse3[1, 2] = new ComplexDouble(0.0, -1.0/30.0);
			Inverse3[3, 1] = Inverse3[1, 3] = new ComplexDouble(-2.0/15.0, 0.0);
			Inverse3[2, 2] = new ComplexDouble(7.0/15, 0.0);

			//            +-     -+
			//            |   1   |
			//            |       |
			//            |  2 I  |
			//            |       |
			// X3 =       |   3   |
			//            |       |
			//            |  4 I  |
			//            |       |
			//            |   5   |
			//            +-     -+
			
			X3 = new ComplexDoubleVector(5);
			X3[0] = new ComplexDouble(1.0, 0.0);
			X3[1] = new ComplexDouble(0.0, 2.0);
			X3[2] = new ComplexDouble(3.0, 0.0);
			X3[3] = new ComplexDouble(0.0, 4.0);
			X3[4] = new ComplexDouble(5.0, 0.0);

			//            +-      -+
			//            |    3   |
			//            |        |
			//            |  48 I  |
			//            |        |
			// Y3 =       |    9   |
			//            |        |
			//            |  60 I  |
			//            |        |
			//            |   15   |
			//            +-      -+

			Y3 = new ComplexDoubleVector(5);
			Y3[0] = new ComplexDouble(3.0, 0.0);
			Y3[1] = new ComplexDouble(0.0, 48.0);
			Y3[2] = new ComplexDouble(9.0, 0.0);
			Y3[3] = new ComplexDouble(0.0, 60.0);
			Y3[4] = new ComplexDouble(15.0, 0.0);
			
			//            +-         -+
			//            |  -6/17 I  |
			//            |           |
			//            |   -16/17  |
			// A3         |           |
			//            |   4/17 I  |
			//            |           |
			//            |    7/17   |
			//            +-         -+

			A3 = new ComplexDoubleVector(4);
			A3[0] = new ComplexDouble(0.0, -6.0/17.0);
			A3[1] = new ComplexDouble(-16.0/17.0, 0.0);
			A3[2] = new ComplexDouble(0.0, 4.0/17.0);
			A3[3] = new ComplexDouble(7.0/17.0, 0.0);
			
		}
		
		#endregion Text Fixture Setup
		
		#region ConstructorNullParameter Test
		
		// Test constructor with a null parameter
		[Test]
		[ExpectedException(typeof(System.ArgumentNullException))]
		public void ConstructorNullParameter()
		{
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(null);
		}
		
		#endregion ConstructorNullParameter Test

		#region ConstructorZeroLengthVector
		
		// Test constructor with a zero length vector parameter
		[Test]
		[ExpectedException(typeof(System.RankException))]
		public void ConstructorZeroLengthVector()
		{
			ComplexDoubleVector dv = new ComplexDoubleVector(1, ComplexDouble.Zero);
			dv.RemoveAt(0);
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(dv);
		}
		
		#endregion ConstructorZeroLengthVector
		
		#region GetVector Test
		
		[Test]
		public void GetVector()
		{
			ComplexDoubleVector dv = new ComplexDoubleVector(new ComplexDouble[]{ComplexDouble.One, 2.0 * ComplexDouble.One, 3.0 * ComplexDouble.One, 4.0 * ComplexDouble.One});
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(dv);
			ComplexDoubleVector dslv = dsl.GetVector();
			Assert.IsTrue(dslv.Equals(dv));
		}
		
		#endregion VectorProperty Test

		#region GetMatrix Test
		
		[Test]
		public void GetMatrix()
		{
			ComplexDoubleMatrix dm = new ComplexDoubleMatrix(3);
			dm[0, 0] = dm[1, 1] = dm[2, 2] = ComplexDouble.One;
			dm[0, 1] = dm[1, 0] = dm[1, 2] = dm[2, 1] = 2.0 * ComplexDouble.One;
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(dm.GetRow(0));
			ComplexDoubleMatrix dsldm = dsl.GetMatrix();
			Assert.IsTrue(dsldm.Equals(dm));
		}
		
		#endregion GetMatrix Test

		#region Decomposition Test 1

		// test the UDL factorisation for case 1
		[Test]
		public void Decomposition1()
		{
			ComplexDoubleSymmetricLevinson dsl = new ComplexDoubleSymmetricLevinson(T1.GetRow(0));

			// get the factors
			ComplexDoubleMatrix U = dsl.GetUpper();
			ComplexDoubleMatrix D = dsl.GetDiagonal();
			ComplexDoubleMatrix L = dsl.GetLower();

			// calculate inverse by the product UDL
			ComplexDoubleMatrix Inv = U * D * L;

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
            	for (int j = 0; i < 3; i++)
				{
                	double E = ComplexMath.Absolute((Inv[i, j] - Inverse1[i, j]) / Inverse1[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			//nunit for mono won't compile this: Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = {0}", MaxError.ToString("E3"));
			Assert.IsTrue(MaxError < 1.0E-15);
		}
		
		#endregion Decomposition Test 1

		#region Decomposition Test 2
		
		// test the UDL factorisation for case 2
		[Test]
		public void Decomposition2()
		{
			ComplexDoubleSymmetricLevinson dsl = new ComplexDoubleSymmetricLevinson(T2.GetRow(0));

			// get the factors
			ComplexDoubleMatrix U = dsl.GetUpper();
			ComplexDoubleMatrix D = dsl.GetDiagonal();
			ComplexDoubleMatrix L = dsl.GetLower();

			// calculate inverse by the product UDL
			ComplexDoubleMatrix Inv = U * D * L;

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
            	for (int j = 0; i < 4; i++)
				{
                	double E = ComplexMath.Absolute((Inv[i, j] - Inverse2[i, j]) / Inverse2[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			//nunit for mono won't compile this: Assert.IsTrue(MaxError < 1.0E-13, "Maximum Error = {0}", MaxError.ToString("E3"));
			Assert.IsTrue(MaxError < 1.0E-13);
		}
		
		#endregion Decomposition Test 2

		#region Decomposition Test 3
		
		// test the UDL factorisation for case 2
		[Test]
		public void Decomposition3()
		{
			ComplexDoubleSymmetricLevinson dsl = new ComplexDoubleSymmetricLevinson(T3.GetRow(0));

			// get the factors
			ComplexDoubleMatrix U = dsl.GetUpper();
			ComplexDoubleMatrix D = dsl.GetDiagonal();
			ComplexDoubleMatrix L = dsl.GetLower();

			// calculate inverse by the product UDL
			ComplexDoubleMatrix Inv = U * D * L;

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 5; i++)
            {
            	for (int j = 0; i < 5; i++)
				{
                	double E = ComplexMath.Absolute((Inv[i, j] - Inverse3[i, j]) / Inverse3[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.5E-15, "Maximum Error = "+MaxError.ToString());
		}
		
		
		#endregion Decomposition Test 3
		
		#region Decomposition Test 4

		// test the UDL factorisation with a singular matrix
		[Test]
		public void Decomposition4()
		{
			ComplexDoubleVector T = new ComplexDoubleVector(4);
			T[0] = T[1] = T[2] = T[3] = new ComplexDouble(1.0, 1.0);
			
			// calculate inverse by the product UDL
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(T);
			dsl.Compute();
			
			Assert.IsTrue(dsl.IsSingular);
		}
		
		#endregion Decomposition Test 4

		#region Determinant Test 1

		// Test the Determinant for case 1
		[Test]
		public void Determinant1()
		{
			// calculate determinant from diagonal
			ComplexDoubleSymmetricLevinson dsl = new ComplexDoubleSymmetricLevinson(T1.GetRow(0));
			ComplexDouble det = dsl.GetDeterminant();
			
			// check results match
			Assert.IsTrue(ComplexMath.Absolute( (det - Det1)/Det1 ) < 1.0E-15);
		}
		
		#endregion Determinant Test 1
		
		#region Determinant Test 2

		// Test the Determinant for case 2
		[Test]
		public void Determinant2()
		{
			// calculate determinant from diagonal
			ComplexDoubleSymmetricLevinson dsl = new ComplexDoubleSymmetricLevinson(T2.GetRow(0));
			ComplexDouble det = dsl.GetDeterminant();
			
			// check results match
			Assert.IsTrue(ComplexMath.Absolute( (det - Det2)/Det2 ) < 1.0E-15);
		}
		
		#endregion Determinant Test 2
		
		#region Determinant Test 3

		// Test the Determinant for case 3
		[Test]
		public void Determinant3()
		{
			// calculate determinant from diagonal
			ComplexDoubleSymmetricLevinson dsl = new ComplexDoubleSymmetricLevinson(T3.GetRow(0));
			ComplexDouble det = dsl.GetDeterminant();
			
			// check results match
			Assert.IsTrue(ComplexMath.Absolute( (det - Det3)/Det3 ) < 1.5E-15);
		}
		
		#endregion Determinant Test 3

		#region SolveVector 1
		
		// solve the linear system for case 1
		[Test]
		public void SolveVector1()
		{
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(T1.GetRow(0));
			ComplexDoubleVector X = dsl.Solve(Y1);
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
				double E = ComplexMath.Absolute((X[i]-X1[i])/X1[i]);
				if (E > MaxError)
				{
					MaxError = E;
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+MaxError.ToString());
		}
		
		#endregion SolveVector 1
		
		#region SolveVector 2
		
		// solve the linear system for case 2
		[Test]
		public void SolveVector2()
		{
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(T2.GetRow(0));
			ComplexDoubleVector X = dsl.Solve(Y2);
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
				double E = ComplexMath.Absolute((X[i]-X2[i])/X2[i]);
				if (E > MaxError)
				{
					MaxError = E;
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-14, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion SolveVector 2
		
		#region SolveVector 3
		
		// solve the linear system for case 3
		[Test]
		public void SolveVector3()
		{
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(T3.GetRow(0));
			ComplexDoubleVector X = dsl.Solve(Y3);
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 5; i++)
            {
				double E = ComplexMath.Absolute((X[i]-X3[i])/X3[i]);
				if (E > MaxError)
				{
					MaxError = E;
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-14, "Maximum Error = "+MaxError.ToString());
		}
		
		#endregion SolveVector 2

		#region SolveMatrix 1

		// calculate inverse by solving linear equations with identity RHS - case 1
		[Test]
		public void SolveMatrix1()
		{
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(T1.GetRow(0));

			// calculate inverse by solving the linear system
			ComplexDoubleMatrix Inverse = dsl.Solve(ComplexDoubleMatrix.CreateIdentity(3));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
            	for (int j = 0; i < 3; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse1[i, j]) / Inverse1[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+MaxError.ToString());
		}
		
		#endregion SolveMatrix 1
		
		#region SolveMatrix 2

		// calculate inverse by solving linear equations with identity RHS - case 2
		[Test]
		public void SolveMatrix2()
		{
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(T2.GetRow(0));

			// calculate inverse by solving the linear system
			ComplexDoubleMatrix Inverse = dsl.Solve(ComplexDoubleMatrix.CreateIdentity(4));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
            	for (int j = 0; i < 4; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse2[i, j]) / Inverse2[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-14, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion SolveMatrix 2
		
		#region SolveMatrix 3

		// calculate inverse by solving linear equations with identity RHS - case 3
		[Test]
		public void SolveMatrix3()
		{
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(T3.GetRow(0));

			// calculate inverse by solving the linear system
			ComplexDoubleMatrix Inverse = dsl.Solve(ComplexDoubleMatrix.CreateIdentity(5));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 5; i++)
            {
            	for (int j = 0; i < 5; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse3[i, j]) / Inverse3[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-14, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion SolveMatrix 3
		
		#region GetInverse 1

		// calculate inverse using GetInverse member - case 1
		[Test]
		public void GetInverse1()
		{
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(T1.GetRow(0));

			// calculate the inverse using Trench's algorithm
			ComplexDoubleMatrix Inverse = dsl.GetInverse();

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
            	for (int j = 0; i < 3; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse1[i, j]) / Inverse1[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion GetInverse 1

		#region GetInverse 2

		// calculate inverse using GetInverse member - case 2
		[Test]
		public void GetInverse2()
		{
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(T2.GetRow(0));

			// calculate the inverse using Trench's algorithm
			ComplexDoubleMatrix Inverse = dsl.GetInverse();

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
            	for (int j = 0; i < 4; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse2[i, j]) / Inverse2[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion GetInverse 2

		#region GetInverse 3

		// calculate inverse using GetInverse member - case 3
		[Test]
		public void GetInverse3()
		{
			ComplexDoubleSymmetricLevinson	dsl = new ComplexDoubleSymmetricLevinson(T3.GetRow(0));

			// calculate the inverse using Trench's algorithm
			ComplexDoubleMatrix Inverse = dsl.GetInverse();

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 5; i++)
            {
            	for (int j = 0; i < 5; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse3[i, j]) / Inverse3[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion GetInverse 3
		

		#region StaticSolveVectorNullParameter 1 Test
		
		// test null parameter
		[Test]
		[ExpectedException(typeof(ArgumentNullException))]
		public void StaticSolveVectorNullParameter1()
		{
			ComplexDoubleVector X = ComplexDoubleSymmetricLevinson.Solve(null, Y1 );
		}

		#endregion StaticSolveVectorNullParameter 1 Test
		
		#region StaticSolveVectorNullParameter 2 Test
		
		// test null parameter
		[Test]
		[ExpectedException(typeof(ArgumentNullException))]
		public void StaticSolveVectorNullParameter2()
		{
			ComplexDoubleVector X = ComplexDoubleSymmetricLevinson.Solve(T1.GetRow(0), (ComplexDoubleVector) null);
		}

		#endregion StaticSolveVectorNullParameter 2 Test
		
		#region StaticSolveVectorDimension Test

		// test mismatching dimensions
		[Test]
		[ExpectedException(typeof(RankException))]
		public void StaticSolveMatrixDimension()
		{
			ComplexDoubleVector X = ComplexDoubleSymmetricLevinson.Solve(T1.GetRow(0), Y2);
		}
		
		#endregion StaticSolveVectorDimension Test
		
		#region StaticSolveVectorSingular
		
		// test with Toeplitz matrix which has a singular principal sub-matrix
		[Test]
		[ExpectedException(typeof(SingularMatrixException))]
		public void StaticSolveVectorSingular()
		{
			ComplexDoubleVector T = new ComplexDoubleVector(4);
			T[0] = T[1] = T[2] = T[3] = new ComplexDouble(1.0, 1.0);
			ComplexDoubleVector X = ComplexDoubleSymmetricLevinson.Solve(T, Y2);
		}
		
		#endregion StaticSolveVectorSingular
		
		#region StaticSolveVector 1 Test
		
		[Test]
		public void StaticSolveVector1()
		{
			ComplexDoubleVector X = ComplexDoubleSymmetricLevinson.Solve(T1.GetRow(0), Y1);
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
				double E = ComplexMath.Absolute((X[i]-X1[i])/X1[i]);
				if (E > MaxError)
				{
					MaxError = E;
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion StaticSolveVector 1 Test
		
		#region StaticSolveVector 2 Test
		
		[Test]
		public void StaticSolveVector2()
		{
			ComplexDoubleVector X = ComplexDoubleSymmetricLevinson.Solve(T2.GetRow(0), Y2);
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
				double E = ComplexMath.Absolute((X[i]-X2[i])/X2[i]);
				if (E > MaxError)
				{
					MaxError = E;
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-14, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion StaticSolveVector 2 Test
		
		#region StaticSolveVector 3 Test
		
		[Test]
		public void StaticSolveVector3()
		{
			ComplexDoubleVector X = ComplexDoubleSymmetricLevinson.Solve(T3.GetRow(0), Y3);
			
			// determine the maximum error
            double MaxError = 0.0;
            for (int i = 0; i < 5; i++)
            {
				double E = ComplexMath.Absolute((X[i]-X3[i])/X3[i]);
				if (E > MaxError)
				{
					MaxError = E;
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-14, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion StaticSolveVector 3 Test

		#region StaticSolveMatrixNullParameter 1 Test

		// test null parameter
		[Test]
		[ExpectedException(typeof(ArgumentNullException))]
		public void StaticSolveMatrixNullParameter1()
		{
			ComplexDoubleMatrix X = ComplexDoubleSymmetricLevinson.Solve(null, Inverse1);
		}
		
		#endregion StaticSolveMatrixNullParameter 1 Test
		
		#region StaticSolveMatrixNullParameter 2 Test
		
		// test null parameter
		[Test]
		[ExpectedException(typeof(ArgumentNullException))]
		public void StaticSolveMatrixNullParameter2()
		{
			ComplexDoubleMatrix Y = ComplexDoubleSymmetricLevinson.Solve(T1.GetRow(0), (DoubleMatrix) null);
		}
		
		#endregion StaticSolveMatrixNullParameter 2 Test
		
		#region StaticSolveMatrixDimension 1 Test

		// test mismatching dimensions
		[Test]
		[ExpectedException(typeof(RankException))]
		public void StaticSolveMatrixDimension1()
		{
			ComplexDoubleMatrix X = ComplexDoubleSymmetricLevinson.Solve(T1.GetRow(0), Inverse2);
		}
		
		#endregion StaticSolveMatrixDimension 1 Test
		
		#region StaticSolveMatrixSingular

		// test with Toeplitz matrix which has a singular principal sub-matrix
		[Test]
		[ExpectedException(typeof(SingularMatrixException))]
		public void StaticSolveMatrixSingular()
		{
			ComplexDoubleVector T = new ComplexDoubleVector(4);
			T[0] = T[1] = T[2] = T[3] = new ComplexDouble(1.0, 1.0);
			
			ComplexDoubleMatrix X = ComplexDoubleSymmetricLevinson.Solve(T, ComplexDoubleMatrix.CreateIdentity(4));
		}
		
		#endregion StaticSolveMatrixSingular
		
		#region StaticSolveMatrix 1 Test

		// calculate inverse by solving linear system with identity RHS
		[Test]
		public void StaticSolveMatrix1()
		{
			// calculate inverse by solving the linear system
			ComplexDoubleMatrix Inverse = ComplexDoubleSymmetricLevinson.Solve(T1.GetRow(0), ComplexDoubleMatrix.CreateIdentity(3));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
            	for (int j = 0; i < 3; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse1[i, j]) / Inverse1[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion StaticSolveMatrix 1 Test

		#region StaticSolveMatrix 2 Test

		// calculate inverse by solving linear system with identity RHS
		[Test]
		public void StaticSolveMatrix2()
		{
			// calculate inverse by solving the linear system
			ComplexDoubleMatrix Inverse = ComplexDoubleSymmetricLevinson.Solve(T2.GetRow(0), ComplexDoubleMatrix.CreateIdentity(4));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
            	for (int j = 0; i < 4; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse2[i, j]) / Inverse2[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-14, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion StaticSolveMatrix 2 Test

		#region StaticSolveMatrix 3 Test

		// calculate inverse by solving linear system with identity RHS
		[Test]
		public void StaticSolveMatrix3()
		{
			// calculate inverse by solving the linear system
			ComplexDoubleMatrix Inverse = ComplexDoubleSymmetricLevinson.Solve(T3.GetRow(0), ComplexDoubleMatrix.CreateIdentity(5));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 5; i++)
            {
            	for (int j = 0; i < 5; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse3[i, j]) / Inverse3[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check the magnitude of error
			Assert.IsTrue(MaxError < 1.0E-14, "Maximum Error = "+ MaxError.ToString());
		}
		
		#endregion StaticSolveMatrix 3 Test

		#region StaticYuleWalkerNullParameter Test
		
		// test Yule-Walker with a null reference
		[Test]
		[ExpectedException(typeof(ArgumentNullException))]
		public void StaticYuleWalkerNullParameter()
		{
			ComplexDoubleVector Y = ComplexDoubleSymmetricLevinson.YuleWalker(null);
		}
		
		#endregion StaticYuleWalkerNullParameter Test
		
		#region StaticYuleWalkerDimension Test

		// test Yule-Walker with order 1 matrix
		[Test]
		[ExpectedException(typeof(ArgumentOutOfRangeException))]
		public void StaticYuleWalkerDimension()
		{
			ComplexDoubleVector T = new ComplexDoubleVector(1);
			T[0] = new ComplexDouble(1.0);
			ComplexDoubleVector Y = ComplexDoubleSymmetricLevinson.YuleWalker(T);
		}
		
		#endregion StaticYuleWalkerDimension Test
		
		#region StaticYuleWalkerSingular
		
		// test Yule-Walker with matrix with singular principal sub-matrix
		[Test]
		[ExpectedException(typeof(SingularMatrixException))]
		public void StaticYuleWalkerSingular()
		{
			ComplexDoubleVector T = new ComplexDoubleVector(4);
			T[0] = T[1] = T[2] = T[3] = new ComplexDouble(1.0, 1.0);
			
			ComplexDoubleVector Y = ComplexDoubleSymmetricLevinson.YuleWalker(T);
		}
		
		#endregion StaticYuleWalkerSingular

		#region StaticYuleWalker 1 Test
		
		[Test]
		public void StaticYuleWalker1()
		{
			ComplexDoubleVector A = ComplexDoubleSymmetricLevinson.YuleWalker(T1.GetRow(0));

			// find the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 2; i++)
            {
				double E = ComplexMath.Absolute((A[i]-A1[i])/A1[i]);
				if (E > MaxError)
				{
                	MaxError = E;
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}

		#endregion StaticYuleWalker 1 Test

		#region StaticYuleWalker 2 Test
		
		[Test]
		public void StaticYuleWalker2()
		{
			ComplexDoubleVector A = ComplexDoubleSymmetricLevinson.YuleWalker(T2.GetRow(0));

			// find the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
				double E = ComplexMath.Absolute((A[i]-A2[i])/A2[i]);
				if (E > MaxError)
				{
                	MaxError = E;
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}

		#endregion StaticYuleWalker 2 Test

		#region StaticYuleWalker 3 Test
		
		[Test]
		public void StaticYuleWalker3()
		{
			ComplexDoubleVector A = ComplexDoubleSymmetricLevinson.YuleWalker(T3.GetRow(0));

			// find the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
				double E = ComplexMath.Absolute((A[i]-A3[i])/A3[i]);
				if (E > MaxError)
				{
                	MaxError = E;
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}

		#endregion StaticYuleWalker 3 Test

		#region StaticInverseNullParameter Test

		// ComplexDoubleSymmetricLevinson.Inverse(ComplexDoubleVector T) with null parameter.
		[Test]
		[ExpectedException(typeof(System.ArgumentNullException))]
		public void StaticInverseNullParameter()
		{
			ComplexDoubleMatrix Y = ComplexDoubleSymmetricLevinson.Inverse(null);
		}

		#endregion StaticInverseNullParameter Test

		#region StaticInverseDimension Test

		// ComplexDoubleSymmetricLevinson.Inverse(ComplexDoubleVector T) with small dimension
		[Test]
		[ExpectedException(typeof(System.RankException))]
		public void StaticInverseDimension()
		{
			ComplexDoubleVector T = new ComplexDoubleVector(1);
			T[0] = new ComplexDouble(1.0, 1.0);
			
			ComplexDoubleMatrix Y = ComplexDoubleSymmetricLevinson.Inverse(T);
		}

		#endregion StaticInverseDimension Test

		#region StaticInverseSingular

		// ComplexDoubleSymmetricLevinson.Inverse(ComplexDoubleVector T) with an ill-conditioned system
		[Test]
		[ExpectedException(typeof(SingularMatrixException))]
		public void StaticInverseSingular()
		{
			// setup an ill-conditioned system (second order principal submatrix is singular)
			ComplexDoubleVector T = new ComplexDoubleVector(4);
			T[0] = T[1] = T[2] = T[3] = new ComplexDouble(1.0, 1.0);
			
			ComplexDoubleMatrix Y = ComplexDoubleSymmetricLevinson.Inverse(T);
		}

		#endregion StaticInverseSingular

		#region StaticInverse1 Test

		// ComplexDoubleSymmetricLevinson.Inverse(ComplexDoubleVector T)
		[Test]
		public void StaticInverse1()
		{
			// calculate the inverse
			ComplexDoubleMatrix Inverse = ComplexDoubleSymmetricLevinson.Inverse(T1.GetRow(0));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 3; i++)
            {
            	for (int j = 0; i < 3; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse1[i, j]) / Inverse1[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}

		#endregion StaticInverse1 Test

		#region StaticInverse2 Test

		// ComplexDoubleSymmetricLevinson.Inverse(ComplexDoubleVector T)
		[Test]
		public void StaticInverse2()
		{
			// calculate the inverse
			ComplexDoubleMatrix Inverse = ComplexDoubleSymmetricLevinson.Inverse(T2.GetRow(0));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 4; i++)
            {
            	for (int j = 0; i < 4; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse2[i, j]) / Inverse2[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-15, "Maximum Error = "+ MaxError.ToString());
		}

		#endregion StaticInverse2 Test
		
		#region StaticInverse3 Test

		// ComplexDoubleSymmetricLevinson.Inverse(ComplexDoubleVector T)
		[Test]
		public void StaticInverse3()
		{
			// calculate the inverse
			ComplexDoubleMatrix Inverse = ComplexDoubleSymmetricLevinson.Inverse(T3.GetRow(0));

			// determine the maximum relative error
            double MaxError = 0.0;
            for (int i = 0; i < 5; i++)
            {
            	for (int j = 0; i < 5; i++)
				{
                	double E = ComplexMath.Absolute((Inverse[i, j] - Inverse3[i, j]) / Inverse3[i, j]);
                    if (E > MaxError)
                    {
                        MaxError = E;
                    }
				}
            }

			// check if test was successful
			Assert.IsTrue(MaxError < 1.0E-14, "Maximum Error = "+ MaxError.ToString());
		}

		#endregion StaticInverse2 Test

	}
	
}
